package fem2.level_set;

import inf.jlinalg.BLAM;
import inf.jlinalg.IVector;
import inf.math.LinAlg;

import java.util.ArrayList;
import java.util.HashSet;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;
import java.util.TreeSet;

import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.LinAlgUtilities;
import math2.MathUtilities;
import fem2.BasisManager;
import fem2.Debugger;
import fem2.Element;
import fem2.IntegrationManager;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshPart;
import fem2.Model;
import fem2.Node;
import fem2.Topology;
import fem2.element.StructuralElement;
import fem2.enu.PropertyName;

/**
 * a crack representation by intersection of two level set
 * 
 * @author hbui
 * 
 */
public class Crack extends Discontinuity {

	private double[] tip;
	private double[] dir;

	private Set<Element> splitElements = new HashSet<Element>();
	private Set<Element> crackTipElements = new HashSet<Element>();

	private double jRadii = 0.0;
	private double propagationStep = 0.0;

	/**
	 * crack constructor
	 * 
	 * @param m
	 *            the associated model
	 * @param tip_start
	 *            start of the tip
	 * @param tip_delta
	 *            displacement vector from start to end of the tip
	 */
	public Crack(Model m, double[] tip_start, double[] tip_delta) {
		super();

		addLevelSet(new LevelSet("psi", m));
		addLevelSet(new LevelSet("phi", m));

		int dim = tip_delta.length;

		this.tip = tip_start.clone();
		LinAlg.add(dim, 1.0, tip_delta, 1.0, tip);

		this.dir = new double[dim];
		double norm = LinAlg.norm2(tip_delta.length, tip_delta);
		LinAlg.scale(dim, 1.0 / norm, tip_delta, dir);

		initialize();
	}

	/**
	 * set the radius to compute interaction integral
	 * 
	 * @param r
	 */
	public void setJIntegralRadius(double r) {
		jRadii = r;
	}

	/**
	 * set step length for crack propagation
	 * 
	 * @param r
	 */
	public void setPropagationStep(double r) {
		propagationStep = r;
	}

	/**
	 * get the step length for propagation
	 * 
	 * @return
	 */
	public double getPropagationStep() {
		return propagationStep;
	}

	/**
	 * get the split elements cut by crack
	 * 
	 * @return
	 */
	public Set<Element> getSplitElements() {
		return splitElements;
	}

	/**
	 * get crack tip element
	 * 
	 * @return
	 */
	public Set<Element> getCrackTipElements() {
		return crackTipElements;
	}

	/**
	 * get crack tip
	 * 
	 * @return
	 */
	public double[] getCrackTip() {
		return tip;
	}

	/**
	 * initialize the crack
	 */
	public void initialize() {
		LevelSet psi = getLevelSet("psi");
		LevelSet phi = getLevelSet("phi");
		Model m = phi.getModel();
		Mesh mesh = m.getMesh();

		/*
		 * update value of level set
		 */
		double a = -dir[1];
		double b = dir[0];
		double c = -(a * tip[0] + b * tip[1]);
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] p = n.getPosition();
			psi.update(n, a * p[0] + b * p[1] + c);
			phi.update(n, (p[0] - tip[0]) * dir[0] + (p[1] - tip[1]) * dir[1]);
		}

		update(psi, phi, m);
	}

	/**
	 * propagate this level set with a distance r and relative angle theta
	 * 
	 * @param r
	 * @param theta
	 */
	public void propagate(double r, double theta) {
		LevelSet psi = getLevelSet("psi");
		LevelSet phi = getLevelSet("phi");
		Model m = phi.getModel();
		Mesh mesh = m.getMesh();

		/*
		 * determine new direction
		 */
		double alpha = computeTipAngle();
		dir[0] = Math.cos(alpha + theta);
		dir[1] = Math.sin(alpha + theta);

		/*
		 * update value of level set
		 */
		double a = -dir[1];
		double b = dir[0];
		double c = -(a * tip[0] + b * tip[1]);
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] p = n.getPosition();
			double phii = (p[0] - tip[0]) * dir[0] + (p[1] - tip[1]) * dir[1];
			// if (phii >= 0) {
			if (phii > 0) {
				psi.update(n, a * p[0] + b * p[1] + c);
			}
			phi.update(n, phii - r);
		}

		/*
		 * determine new tip
		 */
		tip[0] += r * dir[0];
		tip[1] += r * dir[1];

		/*
		 * update split elements and crack tip element
		 */
		update(psi, phi, m);

	}

	/**
	 * @param psi
	 * @param phi
	 * @param m
	 * @throws Error
	 */
	private void update(LevelSet psi, LevelSet phi, Model m) throws Error {
		/*
		 * update the zero level element
		 */
		for (int i = 0; i < m.countElements(); i++) {
			phi.checkZeroElement(m.getElement(i));
			psi.checkZeroElement(m.getElement(i));
		}

		/*
		 * extract split elements and crack tip elements
		 */
		splitElements.clear();
		crackTipElements.clear();

		Set<Element> psiElements = psi.getZeroLevelElements();
		Set<Element> possibleTipElements = new HashSet<Element>();

		for (Element e : psiElements) {
			Node[] nodes = e.getMeshPart().getNodes();

			/*
			 * check if element is splitted
			 */
			boolean isCut = true;
			double psi_min = Double.POSITIVE_INFINITY, psi_max = Double.NEGATIVE_INFINITY;
			for (int k = 0; k < nodes.length; k++) {
				isCut &= (phi.valueAt(nodes[k]) < 0);
				double tmp = psi.valueAt(nodes[k]);
				if (tmp <= psi_min) {
					psi_min = tmp;
				}
				if (tmp >= psi_max) {
					psi_max = tmp;
				}
			}
			isCut &= (psi_min * psi_max <= 0);

			/*
			 * check if element contains crack tip
			 */
			boolean isTip = true;
			double phi_min = Double.POSITIVE_INFINITY, phi_max = Double.NEGATIVE_INFINITY;
			for (int k = 0; k < nodes.length; k++) {
				double tmp = phi.valueAt(nodes[k]);
				if (tmp <= phi_min) {
					phi_min = tmp;
				}
				if (tmp >= phi_max) {
					phi_max = tmp;
				}
			}
			isTip &= ((phi_min * phi_max <= 0) && (psi_min * psi_max <= 0));

			if (isCut & !isTip) {
				splitElements.add(e);
			}

			if (isTip) {
				possibleTipElements.add(e);
			}
		}

		/*
		 * extract true element contains crack tip
		 */
		for (Element e : possibleTipElements) {
			// double[] xi = computeLevelSetIntersection(e);
			// MeshPart mp = e.getMeshPart();
			// if (IntegrationManager.getInstance().checkBoundary(mp.getType(),
			// xi)) {
			// crackTipElements.add(e);
			// }

			if (e.getMeshPart().checkInside(tip)) {
				crackTipElements.add(e);
			}
		}

		if (crackTipElements.size() == 0) {
			throw new Error("could not find crack tip element");
		}

		// Debugger.watch("psi level");
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// Debugger.watch(n.getId() + ": " + psi.valueAt(n));
		// }
		// Debugger.watch("phi level");
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// Debugger.watch(n.getId() + ": " + phi.valueAt(n));
		// }

		// Debugger.watch("tip = ", tip);
		//
		// Debugger.watch("split elements");
		// for (Element e : splitElements) {
		// Debugger.watch(e.getId());
		// }
		//
		// Debugger.watch("crack tip elements");
		// for (Element e : crackTipElements) {
		// Debugger.watch(e.getId());
		// }
	}

	// /**
	// * compute the crack tip by intersecting two level set
	// *
	// * @return
	// */
	// private double[] computeCrackTip() {
	// LevelSet psi = getPrimaryLevelSet();
	// LevelSet phi = getSecondaryLevelSet();
	// Set<Element> psiElements = psi.getZeroLevelElements();
	// Set<Element> phiElements = phi.getZeroLevelElements();
	//
	// /*
	// * find intersection of two sets
	// */
	// Set<Element> intersectedElements = new HashSet<Element>(psiElements);
	// intersectedElements.retainAll(phiElements);
	//
	// /*
	// * search the tip in all intersected elements
	// */
	// for (Element e : intersectedElements) {
	// double[] xi = computeLevelSetIntersection(e);
	// MeshPart mp = e.getMeshPart();
	// if (IntegrationManager.getInstance().checkBoundary(mp.getType(), xi)) {
	// return mp.XAt(xi);
	// }
	// }
	//
	// throw new
	// Error("could not find crack tip by intersecting two level sets");
	// }

	/**
	 * compute intersection of two level sets inside an element
	 * 
	 * @param psi
	 * @param phi
	 * @param e
	 * @return
	 */
	public double[] computeLevelSetIntersection(Element e) {
		LevelSet psi = getLevelSet("psi");
		LevelSet phi = getLevelSet("phi");
		/*
		 * solve for intersection of two level set on crack tip element using
		 * Newton-Raphson procedure
		 */
		MeshPart mp = e.getMeshPart();
		int dim = mp.getDimension();
		double[] xi = new double[dim];
		double[] f = new double[dim];
		f[0] = psi.valueAt(mp, xi);
		f[1] = phi.valueAt(mp, xi);
		double norm_f = Math.sqrt(f[0] * f[0] + f[1] * f[1]);
		int cnt = 0;
		while (norm_f > MathUtilities.tol6 && cnt++ < 50) {
			double[] A = new double[dim * dim];
			double[] tmp1 = psi.gradientAt(mp, xi);
			double[] tmp2 = phi.gradientAt(mp, xi);
			for (int i = 0; i < dim; i++) {
				A[i] = tmp1[i];
				A[dim + i] = tmp2[i];
			}
			LinAlgUtilities.backSolve(BLAM.NO_TRANSPOSE, A, f);
			xi[0] -= f[0];
			xi[1] -= f[1];
			f[0] = psi.valueAt(mp, xi);
			f[1] = phi.valueAt(mp, xi);
			norm_f = Math.sqrt(f[0] * f[0] + f[1] * f[1]);
		}

		if (cnt >= 50) {
			throw new Error("iteration does not converge at element " + e.getId());
		}

		return xi;
	}

	/**
	 * get the crack tip angle
	 * 
	 * @return
	 */
	public double computeTipAngle() {
		return Math.atan2(dir[1], dir[0]);
	}

	/**
	 * compute intersection of tangent level set with element
	 * 
	 * @param e
	 * @return list of intersection points in reference coordinate
	 */
	public ArrayList<double[]> computeIntersectionAt(Element e) {
		MeshPart mp = e.getMeshPart();
		Node[] nodes = mp.getNodes();
		int dim = mp.getDimension();
		LevelSet psi = getLevelSet("psi");

		ArrayList<double[]> points = new ArrayList<double[]>();

		FunctionBasisOnRnToR basis = BasisManager.getInstance().getBasis(mp.getType());
		int[][] edges = MeshPart.getListOfEdges(mp.getType());

		for (int[] ij : edges) {

			Node n1 = nodes[ij[0]];
			Node n2 = nodes[ij[1]];
			double v1 = psi.valueAt(n1);
			double v2 = psi.valueAt(n2);

			if (Math.abs(v1 - v2) < MathUtilities.tol6 && Math.abs(v1) > MathUtilities.tol6) {
				continue;
			}

			double[] xi1 = IntegrationManager.getInstance().getReferenceCoordinate(mp.getType(),
					ij[0]);
			double[] xi2 = IntegrationManager.getInstance().getReferenceCoordinate(mp.getType(),
					ij[1]);
			double[] dxi = xi2.clone();
			LinAlg.add(dim, -1.0, xi1, 1.0, dxi);

			FunctionRnToR basis1 = basis.getBasisFunction(ij[0]);
			FunctionRnToR basis2 = basis.getBasisFunction(ij[1]);

			double t = Math.random();

			double[] xi = new double[] { (1 - t) * xi1[0] + t * xi2[0],
					(1 - t) * xi1[1] + t * xi2[1] };

			double N1 = basis1.valueAt(xi);
			double N2 = basis2.valueAt(xi);
			double f = N1 * v1 + N2 * v2;
			int cnt = 0;
			while (Math.abs(f) > MathUtilities.tol6 && cnt++ < 50) {
				double[] dN1 = basis1.gradientAt(xi);
				double[] dN2 = basis2.gradientAt(xi);

				double df = LinAlg.dot(dim, dxi, dN1) * v1 + LinAlg.dot(dim, dxi, dN2) * v2;
				t += -f / df;

				xi[0] = (1 - t) * xi1[0] + t * xi2[0];
				xi[1] = (1 - t) * xi1[1] + t * xi2[1];
				N1 = basis1.valueAt(xi);
				N2 = basis2.valueAt(xi);
				f = N1 * v1 + N2 * v2;
			}

			if (cnt >= 50) {
				throw new Error("iteration does not converge");
			}

			if (t >= 0.0 && t <= 1.0) {
				points.add(xi);
			}

		}

		/*
		 * remove duplicate entries
		 */
		for (int i = 0; i < points.size(); i++) {
			for (int j = i + 1; j < points.size(); j++) {
				double[] p1 = points.get(i);
				double[] p2 = points.get(j);
				if (LinAlg.diffNorm2(p1.length, p1, p2) < MathUtilities.tol6) {
					points.remove(j);
				}
			}
		}

		return points;
	}

	/**
	 * compute stress intensity factor using interaction integral
	 * 
	 * @return
	 */
	public double[] computeSIFs(Model m) {
		/************************************************
		 * extract all elements lie on the integral path
		 ************************************************/
		Mesh mesh = m.getMesh();

		/*
		 * find nodes lie inside the circle
		 */
		Set<Node> insideNodes = new TreeSet<Node>();
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] p = n.getPosition();

			double r = Math.sqrt(Math.pow(p[0] - tip[0], 2) + Math.pow(p[1] - tip[1], 2));
			if (r <= jRadii) {
				insideNodes.add(n);
			}
		}

		// Debugger.watch("inside nodes");
		// for (Node n : insideNodes) {
		// Debugger.watch(n.getId());
		// }

		/*
		 * if inside nodes contain boundary nodes, then remove it
		 */
		Topology topology = mesh.getTopology();
		if (topology == null) {
			throw new Error("mesh topology is undefined");
		}
		Set<Node> boundaryNodes = topology.getBoundaryNodes();
		insideNodes.removeAll(boundaryNodes);

		/*
		 * find elements lie on the integration path and build the weighting
		 * values
		 */
		Set<Element> pathElements = new TreeSet<Element>();
		Map<Node, Double> weight = new TreeMap<Node, Double>();
		for (Node n : insideNodes) {
			Set<Element> neighbours = topology.getNeighbourElements(n);
			for (Element e : neighbours) {
				boolean isInside = true;
				Node[] nodes = e.getMeshPart().getNodes();
				for (Node en : nodes) {
					isInside &= insideNodes.contains(en);
				}
				if (!isInside) {
					pathElements.add(e);
					for (Node en : nodes) {
						weight.put(en, 0.0);
					}
				}
			}
		}
		for (Node n : insideNodes) {
			weight.put(n, 1.0);
		}

		// Debugger.watch("elements lie on the path");
		// for (Element e : pathElements) {
		// Debugger.watch(e.getId());
		// }
		//
		// Debugger.watch("nodes lie around the path");
		// for (Node n : weight.keySet()) {
		// Debugger.watch(n.getId());
		// Debugger.watch("", weight.get(n));
		// }

		/************************************************
		 * loop through each path element to contribute to interaction integral
		 ************************************************/
		double Imode1 = 0.0;
		double Imode2 = 0.0;

		double alpha = computeTipAngle();
		// CCS = [cos(alpha) sin(alpha); -sin(alpha) cos(alpha)];
		double c = Math.cos(alpha);
		double s = Math.sin(alpha);

		/*
		 * extract necessary constant
		 */
		MaterialModel sampleMaterial = ((StructuralElement) m.getElement(0)).getMaterialModel();
		double mu = sampleMaterial.getMaterialProperty(PropertyName.SHEAR_MODULUS);
		double nu = sampleMaterial.getMaterialProperty(PropertyName.POISSON_RATIO);
		double Eeff = sampleMaterial.getMaterialProperty(PropertyName.YOUNG_MODULUS);
		double km = 0.0;
		switch (sampleMaterial.getState()) {
		case PLANE_STRAIN:
			km = 3 - 4 * nu;
			Eeff /= (1 - nu * nu);
			break;
		case PLANE_STRESS:
			km = (3 - nu) / (1 + nu);
			break;
		default:
			throw new Error("this state is not supported");
		}

		for (Element e : pathElements) {
			MeshPart mp = e.getMeshPart();
			Node[] nodes = mp.getNodes();
			int numNodes = nodes.length;
			int ng = e.getNumIntegrationPoints();
			StructuralElement se = (StructuralElement) e;

			double[] q = new double[numNodes];
			for (int i = 0; i < q.length; i++) {
				q[i] = weight.get(nodes[i]);
			}

			for (int k = 0; k < ng; k++) {
				double[] xi = e.getIntegrationPoint(k);
				double w = e.getIntegrationWeight(k);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);

				/*
				 * compute real stress and convert to crack tip coordinate
				 * system
				 */
				IVector stress = (IVector) se.getStress(k).clone();
				double o11 = c * c * stress.get(0) + s * s * stress.get(1) + 2 * s * c
						* stress.get(2);
				double o22 = s * s * stress.get(0) + c * c * stress.get(1) - 2 * s * c
						* stress.get(2);
				double o12 = s * c * (stress.get(1) - stress.get(0)) + (c * c - s * s)
						* stress.get(2);
				double o21 = o12;

				/*
				 * compute real displacement gradient w.r.t crack tip coordinate
				 * system
				 */
				FunctionRnToR u = mp.getSolution(0);
				FunctionRnToR v = mp.getSolution(1);
				double[] gradu = u.gradientAt(xi);
				double[] gradv = v.gradientAt(xi);
				double u1 = c * c * gradu[0] + s * s * gradv[1] + s * c * (gradu[1] + gradv[0]);
				// double u2 = c * c * gradu[1] - s * s * gradv[0] + s * c *
				// (gradv[1] - gradu[0]);
				double v1 = c * c * gradv[0] - s * s * gradu[1] + s * c * (gradv[1] - gradu[0]);
				// double v2 = c * c * gradv[1] + s * s * gradu[0] - s * c *
				// (gradv[0] + gradu[1]);

				/*
				 * compute gradient of weighting function
				 */
				double[] tmp = new double[2];
				double[] gradq = new double[2];
				double[][] gradN = mp.gGradientsAt(xi);
				for (int i = 0; i < numNodes; i++) {
					for (int j = 0; j < 2; j++) {
						tmp[j] += gradN[i][j] * q[i];
					}
				}
				gradq[0] = c * tmp[0] + s * tmp[1];
				gradq[1] = -s * tmp[0] + c * tmp[1];

				/*
				 * convert integration point to local crack tip coordinate
				 */
				double[] gp = mp.XAt(xi);
				double[] gpLoc = new double[] { c * (gp[0] - tip[0]) + s * (gp[1] - tip[1]),
						-s * (gp[0] - tip[0]) + c * (gp[1] - tip[1]) };
				double r = Math.sqrt(Math.pow(gpLoc[0], 2) + Math.pow(gpLoc[1], 2));
				double theta = Math.atan2(gpLoc[1], gpLoc[0]);

				/*
				 * choose auxiliary state as mode 1 (K1 = 1.0, K2 = 0.0) and
				 * compute interaction intergral
				 */
				double K1 = 1.0;
				double K2 = 0.0;
				double sqr = Math.sqrt(r);
				double ct = Math.cos(theta);
				double st = Math.sin(theta);
				double ct2 = Math.cos(theta / 2);
				double st2 = Math.sin(theta / 2);
				double c3t2 = Math.cos(3 * theta / 2);
				double s3t2 = Math.sin(3 * theta / 2);

				double drdx = ct;
				double drdy = st;
				double dtdx = -st / r;
				double dtdy = ct / r;

				double cAuxStress = Math.sqrt(1 / (2 * Math.PI));
				double cAuxDisp = cAuxStress / (2 * mu);

				// Auxiliary stress component for Mode I loading in x-direction
				double auxStress11 = K1 * cAuxStress / sqr * ct2 * (1 - st2 * s3t2);
				// Auxiliary stress component for Mode I loading in y-direction
				double auxStress22 = K1 * cAuxStress / sqr * ct2 * (1 + st2 * s3t2);
				// Auxiliary shear stress component for Mode I loading
				double auxStress12 = K1 * cAuxStress / sqr * st2 * ct2 * c3t2;
				// Auxiliary shear stress component for Mode I loading
				double auxStress21 = auxStress12;

				// Derivative of auxiliary x-displacement with respect to r
				double du1dr = K1 * cAuxDisp * 0.5 / sqr * ct2 * (km - ct);
				// Derivative of auxiliary x-displacement with respect to theta
				double du1dt = K1 * cAuxDisp * sqr * (-0.5 * st2 * (km - ct) + ct2 * st);
				// Derivative of auxiliary y-displacement with respect to r
				double du2dr = K1 * cAuxDisp * 0.5 / sqr * st2 * (km - ct);
				// Derivative of auxiliary y-displacement with respect to theta
				double du2dt = K1 * cAuxDisp * sqr * (0.5 * ct2 * (km - ct) + st2 * st);

				// Auxiliary displacement gradient of u1 with respect to x
				double auxGradDisp11 = du1dr * drdx + du1dt * dtdx;
				// Auxiliary displacement gradient of u1 with respect to y
				double auxGradDisp12 = du1dr * drdy + du1dt * dtdy;
				// Auxiliary displacement gradient of u2 with respect to x
				double auxGradDisp21 = du2dr * drdx + du2dt * dtdx;
				// Auxiliary displacement gradient of u2 with respect to y
				double auxGradDisp22 = du2dr * drdy + du2dt * dtdy;

				double auxStrain11 = auxGradDisp11;
				double auxStrain12 = (auxGradDisp12 + auxGradDisp21) / 2;
				double auxStrain21 = auxStrain12;
				double auxStrain22 = auxGradDisp22;

				double I1 = (o11 * auxGradDisp11 + o21 * auxGradDisp21) * gradq[0]
						+ (o12 * auxGradDisp11 + o22 * auxGradDisp21) * gradq[1];

				double I2 = (auxStress11 * u1 + auxStress21 * v1) * gradq[0]
						+ (auxStress12 * u1 + auxStress22 * v1) * gradq[1];

				double strainEnergy = o11 * auxStrain11 + o12 * auxStrain12 + o21 * auxStrain21
						+ o22 * auxStrain22;

				Imode1 += (I1 + I2 - strainEnergy * gradq[0]) * w * detJ;

				/*
				 * choose auxiliary state as mode 2 (K1 = 0.0, K2 = 1.0) and
				 * compute interaction intergral
				 */
				// Auxiliary stress component for Mode II loading in x-direction
				K1 = 0.0;
				K2 = 1.0;
				auxStress11 = -K2 * cAuxStress / sqr * st2 * (2 + ct2 * c3t2);
				// Auxiliary stress component for Mode II loading in y-direction
				auxStress22 = K2 * cAuxStress / sqr * st2 * ct2 * c3t2;
				// Auxiliary shear stress component for Mode II loading
				auxStress12 = K2 * cAuxStress / sqr * ct2 * (1 - st2 * s3t2);
				// Auxiliary shear stress component for Mode II loading
				auxStress21 = auxStress12;

				// Derivative of auxiliary x-displacement with respect to r
				du1dr = K2 * cAuxDisp * 0.5 / sqr * st2 * (km + 2 + ct);
				// Derivative of auxiliary x-displacement with respect to theta
				du1dt = K2 * cAuxDisp * sqr * (0.5 * ct2 * (km + 2 + ct) - st2 * st);
				// Derivative of auxiliary y-displacement with respect to r
				du2dr = -K2 * cAuxDisp * 0.5 / sqr * ct2 * (km - 2 + ct);
				// Derivative of auxiliary y-displacement with respect to theta
				du2dt = -K2 * cAuxDisp * sqr * (-0.5 * st2 * (km - 2 + ct) - ct2 * st);

				// Auxiliary displacement gradient of u1 with respect to x
				auxGradDisp11 = du1dr * drdx + du1dt * dtdx;
				// Auxiliary displacement gradient of u1 with respect to y
				auxGradDisp12 = du1dr * drdy + du1dt * dtdy;
				// Auxiliary displacement gradient of u2 with respect to x
				auxGradDisp21 = du2dr * drdx + du2dt * dtdx;
				// Auxiliary displacement gradient of u2 with respect to y
				auxGradDisp22 = du2dr * drdy + du2dt * dtdy;

				auxStrain11 = auxGradDisp11;
				auxStrain12 = (auxGradDisp12 + auxGradDisp21) / 2;
				auxStrain21 = auxStrain12;
				auxStrain22 = auxGradDisp22;

				I1 = (o11 * auxGradDisp11 + o21 * auxGradDisp21) * gradq[0]
						+ (o12 * auxGradDisp11 + o22 * auxGradDisp21) * gradq[1];

				I2 = (auxStress11 * u1 + auxStress21 * v1) * gradq[0]
						+ (auxStress12 * u1 + auxStress22 * v1) * gradq[1];

				strainEnergy = o11 * auxStrain11 + o12 * auxStrain12 + o21 * auxStrain21 + o22
						* auxStrain22;

				Imode2 += (I1 + I2 - strainEnergy * gradq[0]) * w * detJ;
			}
		}

		double[] SIF = new double[] { Imode1 * Eeff / 2, Imode2 * Eeff / 2 };

		Debugger.watch("SIF = ", SIF);

		double K1 = SIF[0];
		double K2 = SIF[1];
		double thetac = 2 * Math.atan(0.25 * (K1 / K2 - Math.signum(K2)
				* Math.sqrt(Math.pow(K1 / K2, 2) + 8)));

		Debugger.watch("theta =" + (thetac * 180 / Math.PI));

		return SIF;
	}
}
